Chapter 7 Beta diversity

load("data/data.Rdata")
load("data/beta.Rdata")

7.1 Captive-bred vs wild

load("data/beta_filtered_capwild.Rdata")
load("data/data_host_filtered.Rdata")

wild_pre_samples <- sample_metadata %>% 
  filter(diet!="Post_grass") %>%
  dplyr::select(sample) %>% pull()

genome_counts<- genome_counts %>% 
  column_to_rownames("genome") %>% 
  select(all_of(wild_pre_samples))%>% 
  rownames_to_column("genome")

sample_metadata <- sample_metadata %>% 
  filter(diet!="Post_grass")

beta_q0n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

genome_counts <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)

beta_q1p <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt %>%
#  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

7.1.1 Permanova analysis

set.seed(1234)

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.004224 0.0042240 7.8055    999  0.003 **
Residuals 22 0.011905 0.0005412                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.008
Wild       0.010583      

adonis2(beta_q0n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.01526961	0.2489907	7.29391	0.001
Residual	22	0.04605641	0.7510093	NA	NA
Total	23	0.06132602	1.0000000	NA	NA

Neutral

betadisper(beta_q1n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.006925 0.0069254 0.5404    999  0.491
Residuals 22 0.281957 0.0128162                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.476
Wild        0.47005      

adonis2(beta_q1n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	1.538238	0.2539205	7.487472	0.001
Residual	22	4.519715	0.7460795	NA	NA
Total	23	6.057953	1.0000000	NA	NA

Phylogenetic

betadisper(beta_q1p$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.007671 0.0076708 1.8217    999  0.179
Residuals 22 0.092637 0.0042108                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.187
Wild        0.19084      

adonis2(beta_q1p$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.3720264	0.3082818	9.804858	0.001
Residual	22	0.8347476	0.6917182	NA	NA
Total	23	1.2067741	1.0000000	NA	NA

Functional

betadisper(beta_q1f$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.003070 0.0030698 1.7361    999  0.171
Residuals 22 0.038901 0.0017682                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.176
Wild        0.20119      

adonis2(beta_q1f$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.01079711	0.187686	5.083124	0.103
Residual	22	0.04673041	0.812314	NA	NA
Total	23	0.05752752	1.000000	NA	NA

7.1.2 Beta diversity plots

7.1.2.1 Richness

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet", "Wild" = "Wild"))+
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(face = "bold", size = 12),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 12),
      legend.title = element_text(size = 14),
      legend.position = "right", legend.box = "vertical"
    ) +labs(color='Origin')

7.1.2.2 Neutral

7.1.2.3 Phylogenetic

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.1.2.4 Functional

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
  geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme_classic() +
  theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.2 Before and after grass

load("data/beta_filtered_pre_post.Rdata")

genome_counts <- genome_counts_filt
post_pre_samples <- sample_metadata %>% 
  filter(diet!="Wild") %>%
  dplyr::select(sample) %>% pull()
genome_counts<- genome_counts %>% 
  column_to_rownames("genome") %>% 
  select(all_of(post_pre_samples))%>% 
  rownames_to_column("genome")
sample_metadata <- sample_metadata %>% 
  filter(diet!="Wild")

beta_q0n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

genome_counts <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)

beta_q1p <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

7.2.1 Permanova analysis

set.seed(1234)

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0000047 0.00000467 0.0151    999  0.888
Residuals 22 0.0067978 0.00030899                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                 0.885
Post_grass   0.90328           

adonis2(beta_q0n$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.0006002558	0.030093	0.6825872	0.785
Residual	22	0.0193464336	0.969907	NA	NA
Total	23	0.0199466894	1.000000	NA	NA

Neutral

betadisper(beta_q1n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.001983 0.0019829 0.1447    999  0.686
Residuals 22 0.301550 0.0137068                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                 0.686
Post_grass   0.70733           

adonis2(beta_q1n$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.147998	0.03074149	0.6977632	0.946
Residual	22	4.666278	0.96925851	NA	NA
Total	23	4.814276	1.00000000	NA	NA

Phylogenetic

betadisper(beta_q1p$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.001772 0.0017715 0.4816    999  0.468
Residuals 22 0.080927 0.0036785                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                 0.465
Post_grass   0.49497           

adonis2(beta_q1p$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.02764442	0.04145548	0.951464	0.59
Residual	22	0.63920164	0.95854452	NA	NA
Total	23	0.66684606	1.00000000	NA	NA

Functional

betadisper(beta_q1f$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.008352 0.0083519 3.722    999  0.037 *
Residuals 22 0.049367 0.0022439                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                  0.03
Post_grass  0.066699           

adonis2(beta_q1f$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.002870065	0.03722968	0.8507252	0.452
Residual	22	0.074220722	0.96277032	NA	NA
Total	23	0.077090787	1.00000000	NA	NA

7.2.2 Beta diversity plots

7.2.2.1 Richness

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.2.2.2 Neutral

7.2.2.3 Phylogenetic

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.2.2.4 Functional

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')